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I. INTRODUCTION 

Calcium ions are ubiquitous mediators of metabolic 
change in many cellular systems. In the cortex, it is 
believed to be the primary chemical signal for the induc- 
tion of bidirectional synaptic plasticity 0,0, H3|, which is 
widely accepted to be the basis for many forms of learn- 
ing, memory and development. 

Experimentally, bi-directional synaptic plasticity can 
be induced by various different induction protocols [|| 
[Tot IT^ I , including the recently characterized spike-time- 
dependent plasticity (STDP) fl |22|. In this induction 
method, if a presynaptic action potential occurs within a 
window of a few tens milliseconds before a postsynaptic 
action potential (At > 0), long-term potentiation (LTP) 
is elicited; if the order is switched (At < 0), long-term 
depression (LTD) happens. Interestingly, STDP and 
many other induction protocols share a common prop- 
erty, which is their dependence on the amount of inte- 
grated activation of the N-methyl-D-aspartate receptors 
(NMDAR) during conditioning. NMDAR are natural co- 
incidence detectors; their activation relies both on the ef- 
ficiency of neurotransmitter release from the presynaptic 
cell, and on the level of depolarization of the postsynap- 
tic cell. Most remarkably, they are permeable to calcium 



ions, which suggests a functional link between calcium 
influx through NMDAR and the induction of LTP and 
LTD. 

Recently, several dynamical models have been pro- 
posed to account for STDP 0, El H E3 , the relation- 
ship between STDP and rate-based models has been ex- 
amined using averaging methods (15L [29j , and there have 
been attempts to account for receptive field formation 
on the basis of STDP |33|. We have recently proposed 
a unified theory of synaptic plasticity [2^ based on the 
following assumptions: (1) the level of intracellular cal- 
cium concentration determines the sign and magnitude of 
synaptic plasticity 0, 0, 0] : when calcium falls below a 
threshold 0d, no plasticity occurs, when it falls between 
0d and Op, 6, p > 0d, LTD is induced, and for calcium 
above 9 p , LTP is induced; (2) the relevant sources of 
calcium are the NMDA channels; and (3) dendritic back- 
propagating action potentials (BPAPs) contributing to 
STDP have a slow "after-depolarizing" tail component. 
We have shown that these simple assumptions are suffi- 
cient to account for the various experimental plasticity- 
induction protocols. In addition, this model has also pro- 
duced previously uncharacterized predictions, such as: 
(a) the shape of the STDP learning curve should be 
frequency-dependent; and (b) there should be a novel 
form of spike time-dependent LTD for At larger than the 
LTP-inducing intervals. Recent experimental results are 
consistent with these predictions |24Ll30|. 
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In this paper, we present a systematic calculation of 
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the calcium transients as a function of pre- and post- 
spike timings, as well as of the neuronal firing rates. In 
section |H] we derive the solution for the mean calcium 
transients, in the simplest case where the BPAP consists 
of a single exponential, and show that its dependence 
on At is rather intuitive. A more realistic approach is 
taken in section II I II where the BPAP is composed of a 
sum of two exponentials with different time constants. 
This more complex assumption is in closer agreement 
with our previous work 123 and is consistent with ex- 
perimental observations |20|. The analysis shows that 
the slow component contributes for different calcium lev- 
els in the baseline condition (pre only) and the post-pre 
condition (At < 0), while the fast component sharpens 
the difference of the calcium levels in the transition be- 
tween post/pre (At < 0) and pre/post (At > 0) condi- 
tions. The sharpness of this transition is shown to be 
further enhanced by the appropriate choice of the cal- 
cium decay time-constant. For completeness, in section 
IIVI we study the effects of calcium-accumulation due 
to repetitive firing which affects the rate-dependence of 
some plasticity-induction protocols including the STDP 
learning window. Finally, in section we take into ac- 
count the stochastic properties of synaptic transmission 
and estimate the trial-by-trial variability of the calcium 
transients. This variability can be significant if the num- 
ber of postsynaptic NMDAR in each spine is small. Our 
analysis demonstrates that this variability depends At; 
increasing with At for At > 0. These results could signif- 
icantly alter the form of the STDP curve under different 
experimental conditions. 



In this paper we consider the NMDA channels as the 
sole source of calcium currents. These channels are volt- 
age dependent and have a long open time that can 
usually be described as a sum of two exponentials 0, 
thus its currents can be generally described as: 



I(t,{t pre },{t post })=gH(V)f, (3) 

where {t pre } and {t post } are the sets of times of pre- and 
postsynaptic action-potentials (APs) that occured prior 
to time t and g is the mean total conductance of the pop- 
ulation of NMDAR in a synapse. The 7i-function repre- 
sents the voltage dependence of calcium influx through 
NMDAR, and / is the fraction of NMDAR in open 
state. It is easy to see that V = V(t, {t pre }, {t post }) 
and f = f(t, {t pre }), where the time dependence of / is 
a double exponential. 

In this section, however, we will focus on the simple 
scenario of a single pair of pre and postsynaptic APs, 
where the dynamics of the NMDAR and the BPAP are 
described as single exponentials with time constants r N 
and r B , respectively. We will leave the more realistic 
two-component-BPAP case to the next section. Unless 
otherwise mentioned, we replace the stochastic variable 
/ by its mean without change of notation. Thus, for a 
single presynaptic spike at time t pre , we can write, 



f(t,t pre ) = nO(t-t pre )e- {t ~ tPr ° )/TN , (4) 



II. ANALYSIS OF CALCIUM INFLUX 
THROUGH NMDAR 

The postsynaptic intraspinal calcium ions are removed 
from these compartments both through diffusion and 
through binding to calcium buffers. Let these phenomena 
be characterized by a single time constant r, and assume 
that the calcium dynamics follow a first-order linear dif- 
ferential equation, 



d ^=I(t)- 1 -[Ca}. 
at t 



(1) 



where [Ca] is the intracellular calcium concentration, I(t) 
is the inward calcium current, and r is the time constant 
of the passive decay process. This simple dynamics is to 
the first approximation consistent with experiments plj 
|26j. Given a initial condition [Ca](0), equation ^ has a 
solution of the form: 



[Ca}(t) 



e^I^dt' +[Co](0) 



(2) 



Notice that if the calcium current is a superposition of 
separate current sources, the calcium concentration can 
also be decomposed into such a sum. 



where O is the Heaviside function and /i is the probabil- 
ity that a channel will open given a presynaptic AP. Let 
us assume that the postsynaptic depolarization is domi- 
nated by the BPAP. Therefore, V in equation [3] can be 
written as: 



V(t,t post ) = V R + BPAP(t) 

= V R + V B Q(t - t post )e-^- fP 



*)/tb 



(5) 
(6) 



where Vr is resting membrane potential and Vb is the 
magnitude of the BPAP due to a postsynaptic AP oc- 
curing at time t post . In general, 7i is a non- linear and 
non-monotonic function of V. For tractability, we will 
make a major simplification by linearizing the 7i-function 
in the interval [Vr, Vb], 



H(V) = a + bV. 



(7) 



The original TL function used and its linear fit are dis- 
played in figure lTUl In appendix^] we assess the validity 
of this approximation. By using equations 01 through 
equation |3| becomes: 
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I(t,t pre ,t post ) 

=g (a + bV R + bV B Q{t - ^*) e -(*-W)/nA 

=W (a + bV R ) 9(f - t Pre^ e -{t-t"n/rN 

+ gfibV B 6(t - t posf )@(t - t pre ) 

x e -(^t P ")/TN e -(t-t P ° at )/TB 



=P re (t) +I +/ -{t). 



(8) 



The current separates into two components, one that 
depends only on the timing of the presynaptic AP (IP re ) 
and a second, associative term, that depends on the 
relative timing between the pre- and postsynaptic APs 
(/ + /~). The associative term will give rise to the differ- 
ences between pre/post and post/pre conditions, which 
is essential to the STDP learning. 

Let At = t post — t pre . It is straightforward to show 
that: 



/+/" = 



7 ^eafc e ( < - tP"«*) e -(*- t ""'Vri } if At > 



tpeak^it ~ t pre )e- (t - tVT °^ r \ otherwise, 



(9) 



where t x 1 = r B 1 + r N \ and I+ ak = /+ afe (At) and 

I~ = I" (At) denote the peak magnitudes of the pre/post 
component of the calcium current for At > and At < 
respectively: 

^(At) = gfibV B e- At ^ (10a) 
I- ak (At)=gf,bV B e At ^. (10b) 

The dependence of the peak current on At is quite 
simple (FIG. QJ. For At > it decays exponentially 
with the time constant of the NMDAR open probability 
r N ; whereas for At < it decays with r B , the width 
of the BPAP. It is straightforward to see that, if At < 
yields calcium levels corresponding to LTD, and if to each 
calcium concentration is associated a degree of plasticity, 
there should be a At > region that also corresponds to 
LTD. The difference is that, since r B << r N , 9d in the 
post /pre region is reached for values of At much closer to 
zero than it is for the pre/post region. Such difference will 
be even greater in the case where the BPAP's dynamic 
is described by two time constants instead of one (see 
section ITTT|) . 

Using equations |S1 and 13 we can rewrite equation |21 in 
terms of each of the components of the calcium current, 



[Ca]{t) = [Ca] pre (t) + Q(At)[Ca] + (t) 

+ e(-At)[Ca]-(t) + [Co](0) 



-t/r 




j f -A t/T 

veak 



50 100 150 200 
A t (ms) 

FIG. 1: Associative component of the peak calcium current 
as a function of At. gb = 1.5 x 10" 3 , fx = 0.8, Vb = 60 mV, 
r N = 100 ms and r B = 20 ms. 



with 



[Ca] pre (t) = e- t/T J I pre (t')e t,/T dt' 

[Ca]-{t) = e~ t/T J I-{t')e t,/T dt'. 
Integrating explicitly, we have, 

[Ca] pre (t) = gfi(a + bV R )r 2 Q{t - t pre ) 



-(t-t pr ")/TN 



lCa] + (t)=I+ ak T 3 e(t~t p ° st ) 



L-{t-t»°°* 



)/n _ e -(t-t post )/r 



[Ca]-(t)=I- k r 3 Q(t-t^) 



(12a) 
(12b) 
(12c) 

CE 

(tHcli 



)/n _ e -{t-t vre )h 



where r 2 1 = r 1 — r N 1 and r 3 1 = r 1 — r N 1 — t b 1 . 
Examples of calcium transients are displayed in FIG.|21 



where [Ca] pre/post = [Ca 



(11) 



pre + [Ca}+ and [Ca] post/pre = 
[Ca] pre + [Ca]" . Naturally, the [Ca] pre term is identical 
for the pre-post (FIG.|2^, At = 10 ms) and the post-pre 
(FIG. Ep, At = -10 ms) conditions (dashed line). The 
difference between these conditions lies in the associative 
term (dash-dot line). [Ca] + raises at the point t = t post 
and [Ca]" raises at t = t pre , their relative amplitudes 
being determined by I^~ k (FIG. Q}. 



III. A BACK-PROPAGATING ACTION 
POTENTIAL WITH TWO COMPONENTS 

The model described in the previous section assumes 
a BPAP with a single exponential falling phase. How- 



4 




FIG. 2: The dynamics of the total calcium (solid line) is composed of two different contributions: the calcium due to the 
pre-spike alone (dashed line) and the calcium due to the association of pre- and postspikes (dash-dot line), ag = 1.03a;10 _ , 
Vr = — 65m V, r — 50 ms, [Co](0) = and the onset of t pre is at time zero. The remaining parameters are as in FIG. a) 
At = 10 ms, b) At = -10 ms. 



ever, there are experimental indications that an addi- 
tional slow after-depolarizing com pon ent exists in some 
cell bodies ^lj an d dendrites [rH |2C| . The dynamics of 
the BPAP defines to what extent the postsynaptic spike- 
timing interacts with the presynaptic one; in particular, 
a slower decaying process ensures that such interaction 
spans for a period of tens of milliseconds, as required by 
the STDP learning rule. In this section, we assume a 
BPAP composed by a sum of two exponentials, with, re- 
spectively, a slow t| and a fast t£ time constants. Equa- 
tion should therefore be re- written as: 



v(t, t post ) =v R + v B e(t - t post ) 



x v J e 



f p -{t-t post )/ri +v s e -(t-t"" t )/r 



where v s , , tJ and r£ are the relative amplitudes and 
time constants of the slow and the fast components of 
the BPAP respectively, (v s ,v^) > 0, v s + = 1. Analo- 
gously to the previous section, the calcium currents from 
equation [5] can be separated into two components, I pre 
and T pre is identical to what we have derived be- 

fore, whereas / + /~ reads: 



/+/"(*) = 



®(t-t post K ak 



,/ p -(t-* post )/r/ 



+v e 



-(t-t"° st )/r= 



if At > 



Q{t-t pre ) 



r-f f -(t-i»">)/r/ 
1 peak u e 



+J- s ,v s e^( t -* P " e )/ r i 

1 peak 



otherwise 

0) 



where Ip eak is as defined in equation llOal I p / ak — 



gfj,bV B e At/T B, I pe s ak = gfj,bV B e At/T B. The calcium cur- 
rent IP re /p° st peaks at t = t post for At > with am- 



plitude Ip eak , and at t = t pre for At < with am- 
plitude r eak = v s I 



peak 



V ^ Ipeak- 



Thus » T peak decays 
more abruptly than its single-component-BPAP counter- 
part (FIG. Eli). 

We can again integrate each component of the calcium 
currents separately. We find for this case: 



[Ca] pre (t) = gfx(a + bV R ) T2 Q(t - t pre ) 



((Hal ) 



x j e -(t-« prc )/T N _ e -(t-t pr ")/r 



[Ca] + (t)=I+ ak @(t-t post ) C20) 
x {v f rl (e-^- tPO,t y-( - e -(*-* P ° st )/-) 
+v s T s i-(t-t*°**)M _ e -(t-f~* 

[Co\-(t) = e(t-t pre ) pel ) 



+r s v s rl (e 



where (r/)" 1 = (ri)" 1 + t~\ (t^ 1 = r" 1 - - 
(tI)^ 1 , and rf and r| are defined similarly to r|, but 
with r| replacing . 

The effects of the BPAP slow component on the cal- 
cium transients is shown in FIG.EJa and0]i. Since both 
[Ca} + and [Ca]~ have the same linear dependence on 
the magnitude v s of this slow component, it primarily 
contributes to the separation between the pre-only con- 
dition (diamond), and the associative conditions (circle 
and square). 

A longer decaying tail for the postsynaptic variable 
provides a longer range of time-interval in which inter- 
action with the presynaptic variables is possible. Thus, 
the dynamics of the presynaptic variable, such as the 
calcium time constant r, will also influence the relative 
magnitudes of the calcium transients for the pre-post and 
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(a) (b) 




A t (ms) time (ms) 

FIG. 3: a) Associative component of the peak calcium current as a function of At, for v f = 0.75, v s = 0.25, t( = 3 ms and 
i~b = 35 ms. The dotted line shows the previous, single BPAP component result, for comparison, b) Total calcium transients 
for At = 10 ms (dashed line), -10 ms (dashed-dot line) and [Ca] pre (solid line). All the other parameters are as in FIG.0|5| 



post-pre conditions. Indeed, the combination of a faster 
calcium kinetics with the different BPAP time courses 
enhances the difference between these magnitudes (FIG. 
0Jd). Therefore, a two-component-BPAP model produces 
a sharper transition of the learning curve between the 
post-pre LTD and the pre-post LTP. Initial estimates of 
the calcium time constant in the spines were of the or- 
der of 70-125ms 0,H^|. However, these values are ob- 
tained through calcium imaging, whose results may de- 
pend on the kinetics of the calcium indicators used, lim- 
iting the accuracy of the estimates. Recent experimental 
techniques, in which the effect of the calcium buffers are 
minimized, have lead to time constants of as low as 20 
ms El. 



IV. RATE-DEPENDENCE OF THE CALCIUM 
TRANSIENTS 

In the previous sections we have analyzed the calcium 
transients evoked by one pair of pre and postsynaptic 
spikes, these have significant implications for STDP at 
low frequency. However, STDP at low frequency is not 
the sole method for inducing synaptic plasticity, and in 
some cases even STDP requires that the pairs of pre 
and postspikes are delivered above a certain frequency 
[22j,|30|. Thus it is important to analyze the contribution 
of spike delivery frequency z/, to the calcium concentra- 
tion, in addition to the relative timing effect. We analyze 
the case in which both pre and postsynaptic neurons fire 
regularly at the same frequency u, the onset of the post- 
synaptic spike-train being shifted from the presynaptic 
spike train by At < 1/v. Notice, however, that in the 
most general case the number of presynaptic and postsy- 
naptic spikes as well as their timing differences will vary. 

Recall that the NMDA current in equation [3] is the 
product between the fraction of opened channels f(t) and 
the linearized Mg-block function H(V{t)) = a + bV(t), 
where V(t) = Vr + BPAP(t) omits depolarization due 



to the excitatory postsynaptic potentials (EPSPs). For 
a sequence {i^ re , tf re , t^ e } of presynaptic spike-times 



post .post 



and {t\ os , t 
can write: 



■,£jJJ s } of postsynaptic spike-times, we 



/(*) 



N 

£ 

71=1 



A, 



exp 



t - iE™ 



e(t-c e ), (13) 



M-l r 



BPAP(t) = V B J2 



exp 



f _ fPost \ 

t —^- e(t-c 5t ) 

T B J 



xe(-t + Ci) 



(14) 



exp 



t - 1 



post 
M 



e (t - tT) 



where N and M are the total number of pre and post- 
spikes before time t, and A n is the increase of the fraction 
of opened channels upon each pre-spike n. Let f(t) be 
at level B n immediately before the pre-spike n occurs. 
If the NMDA opening is proportional to the amount of 
non-opened channels, then A n — /i(l — B n ), where [i is 
the fraction of previously closed channels that open due 
to the presynaptic spike, /j, < 1. Writing B n explicitly 
for each n, it is easy to see that it satisfies the following 
expression: 



n p -n/vTti 



1 - (1 - M)"e 
1 - (l-^)e- 1 A 



and 



Boo = km B„ = 



fie 



-l/i/TN 



1 - (1 - /i)e- 1 /^N 



(15) 



(16) 
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FIG. 4: a) Peak calcium transients as a function of the slow component of the BPAP v B for the pre-before-post (circle), post- 
before-pre (square) and the pre-only (diamond) conditions, b) Relative magnitudes of the peak calcium transients, between 
pre and post-pre (circle), and post-pre and pre-post (square) conditions, as a function of the calcium time-constant r. All 
parameters are as in FIG.|3] 



With this, equation IT3I can be written as: 

t - tP re 



N-l 



/(*) = £] (B n + A n )exp 



71=1 



xe(t-c e )e(-i + C+ e i) 



(A N + B N ) exp - 



t - t T 



Q(t-t 



N 



N-l 



n ^ 

x9(t-ae(-t+Ci) 



t - 



(17) 



+ [B N (1 -fji)+fi] exp - 



We can define, analogously to equation and the ex- 
pressions ^a] through ^c] from section the following 
variables: I pre = g(a + bV R )f(i), = gbBPAP(t)f(t) 
and: 



[Ca] pre (t) = e- f / T J P^it'y'^dt', (18a) 

[Ca] +/ -{t)=er t/T J IP re /P° st (t')e t ' /T dt'. (18b) 
Using equation El equation I18al reads: 

[Ca}P re = ge- l l T (a + bV R ) r 2 

( N-l 

7/--l)£[B n (l- M )+ M ]e^ 

v ' n=l 

+ (e t/T * - e*""/" 2 ) [B N (1 -n) + ft] 



x < e 



(19) 



We use the approximation B n — Boo, and write the 
time indexes with respect to the timing of the first spike: 



re _ .pre 



n-l 



.pre .pre , L^J ^ 

N 1 ~ 



and 



N = [(t - if>J + 1, 



(20a) 

(20b) 
(20c) 



where |_ J , is the floor operation. Evaluating the sum, we 
have: 



[Ca]P re =ffe"'/> + bV R )r 2 [B M (1 - a) + a) 



; i/(- 2 ) _ A e *r/ 



l-exp(^) 



exp 



.pre 
l N 



t 

T-2 



!-ex P (^) 
«xpf^i (21) 



where TV and t p ^ e are defined as in equations I20bl and 
33 



To calculate the contribution to the calcium transients 
due to the interaction of the pre and postsynaptic spikes, 
we substitute the expression 1141 and 1171 into equation 
ILSbl Let C° st = *iT + At and i^ st = t p ^ e + At. 
The indexes m and M will always refer to the post- 
spikes while n and N will always refer to the pre- 
spikes; for simplicity, we drop the superscripts so that 
tf 1 re — t n and tP£ st — t m . Because the depolarization 
due to each post-spike does not build up, the product 
between the sums given by expressions 1171 and 1141 will 
yield only four terms, due to the interaction of the pairs 
{(n,m - 1); (n,m); (N,M - 1); (N,M)}. These terms 
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jf f(t')BPAP{t')dt' = V B - /i) + M] I 



E 

n=l 



exp 



^m — 1 


^ t' 






f' - t 


^ t' 








(22d) 



where £* = min(i, £m) and we have used that B n = Bn = 
Boo. Recalling expressions l20aJ through [2Ucl one can eas- 
ily perform the integrations I22al through I22dl The re- 
maining sums are finite power series, therefore they con- 
verge. Thus, the calcium concentration due to the inter- 
action between pre and postsynaptic spikes will be: 



Plots of these expressions are shown in FIG. As 
expected, the overall calcium is frequency and timing- 
dependent. There is temporal integration of calcium lev- 
els, so that higher frequency stimulations build up the in- 
tracellular calcium concentration (FIG. [5]}). In addition, 
the Ai-dependence shown in previous sections (FIG. EJd) 
for single pairs of pre and post spikes is retained. The 
amount of temporal integration will naturally depend 
on the specific values of r, r N and r B , and so will the 
difference between the pre-post and the post-pre condi- 
tions, as analyzed before. For example, slower dynamics 
would result in more time integration at moderate fre- 
quencies. However, given these constants, it is possible 
to set a threshold of calcium concentration between LTD 
and LTP (or between no-plasticity and LTD) that would 
correspond to low- and high-frequency STDPs, respec- 
tively. 

We now extract an expression for the dependence of 
calcium transients on At and on v in the limit t — > oo. 
Let t = t p ^ e + 5, where the new variable 8 tracks the time 
since the last presynaptic spike, and define 71 = g(a + 
&Vr)t 2 [-Boo(1-m) + m] and 72 = 0&Vbt 3 [.B oo (1 - fi) + fx[. 
For t — > 00 and At < S, we have: 



[Ca}+/- = ge- f / T bV B [B M (1 - fi) + fi] 

If2la1 + [2lb1 + f22"c1 + [22d| (23) 

with each of the terms being: 





[At 






exp 






V T 2 





x e 



exp 



At - l/v 



l-c*p(^) 
1-«*p(£) 




l/l/T2 p l/UT 



e l/vr _ 1 



e l/ur _ I 

e At/r 3 _ 1 



,1/UT _ ^ 



(24) 





(At 






T3 


exp i 


--!-) 













x e 



tfVr 



exp I — — 



exp 



x exp 



T3 



T3 



exp 
x exp 



(*), 

\ n T3 ) 
f At-l/tA 



l-cxp(^) 
. l-exp(^) 

At-l/w s 

fl'l' 



if (t^ re + At) < t 



exp 



if (i^ e + At) > t. 



exp 
x exp 



+ - exp (^) 



At 

TN 



0, 



if (4 re + At) < t 
if (t% e + At) > t. 



And for At > 5: 



[Ca}(6) = lie -^ T \e s/T2 + 



,1/W 2 _ gi/^ 1 " 



e l/»/T _ 1 
; At/T 3 _ j 



-1/i/tb 



o1/vt 3 _ p At/r 3 



,1/vt _ I 



( e V-3_^g-l/^| 



(25) 



These expressions show how the peaks of the calcium 
transients depend on the timing and on the frequency 
of the pre and postsynaptic spikes. As the frequency 
increases (FIG. EJl) the Ai-dependent curves move up 
due to temporal integration, and the difference between 
pre-post and post-pre (FIG.IHb) becomes smaller. 
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FIG. 5: a) Rate-dependence of the calcium transients as derived in expressions I22al through |22dl for At = 10 ms at 5 Hz (solid 
line) and 10 Hz (dashed line), b) Calcium transients at 5 Hz for At = +10 ms (solid line) and At = —10 ms (dotted line). All 
parameters as in FIG.0 




FIG. 6: a) Peak calcium transients as a function of the timing between pre and postsynaptic spikes for three different 
frequencies: 1 Hz (circle), 5 Hz (square) and 10 Hz (triangle), b) Peak calcium transients as a function of the frequencies for a 
pre-post (square) and a post-pre (circle) conditions. 



V. STOCHASTIC PROPERTIES OF CALCIUM 
INFLUX 

In previous sections we have calculated the mean cal- 
cium transients under various conditions and assump- 
tions. However, synaptic transmission is a random pro- 
cess and the trial-by-trial calcium transients deviate from 
the mean. One source of variability is the stochasticity 
of the presynaptic neurotransmitter release. A failure 
of release at a given synapse would result in eliminating 
the postsynaptic calcium transient at that synapse. Al- 
though this source of variability can be significant, its 
consequences are simple to predict, because the absence 
of release will result in no synaptic plasticity in the spe- 
cific spine. Another source of variablity is the stochas- 
ticity of the postsynaptic opening of NMD A channels. 
Because the number of NMDAR at each synapse can be 
small (« 10) [2f|, these fluctuations could be consider- 
able. 

In this section, we calculate the variability due to Z 
NMDAR in the postsynaptic terminal. We assume the 



following simple Markov model: the state vector P„(i) 
denotes the probability of being in each one of the n 
possible states and the transition matrix R n represents 
transition probabilities from state to state. For a Markov 
processes, the evolution equation has the form: 

= RnPn, (26) 

at 

with the solution: 

P n (i) = P n (0) exp Qf R n dt^j , (27) 

where P n (0) is the initial condition. 

In the simplest model the process has only two states, 
unbound u and open-bound o: 

fciG 

u^o (28) 
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where G is the glutamate concentration, k\ and are 
the forward and backward time constants respectively. 

This model is clearly simplified; typically, a five-state 
model is used to account for NMDAR kinetics. However, 
it is sufficient for describing the single exponential decay 
kinetics we have assumed for the NMDAR since, during 
decay, G = (see equation 1311 below). The transition 
matrix for this model is: 



R.2 = 



-kiG k-i 
kiG -k- x . 



where G is the neurotransmitter concentration. We as- 
sume that G has a constant, non-zero value for only a 
limited amount of time. We are primarily interested in 
the falling phase (G = 0) of the NMDAR current. 

We denote by p u and p° the probability of being in the 
unbound and open states respectively. Since p u + p° = 1 
the differential equation for p° reads: 



p° 



k 1 G-p°{kiG + k^i) 



(30) 



For a constant value of G, the solution of this equation 



is: 



P°(t) 



kiG 



kxG + k-i 



Ce 



-(fciG+fc_i)t 



(31) 



where the constant C is determined by the initial condi- 
tions. 

During the falling phase (G — 0) the time constant 
has the value r N = 1/k—i. The solution takes the form 
p°{t) = /ie~*/ TN , where [i can now be defined more pre- 
cisely as the open probability at the start of the falling 
phase. The total calcium current through Z identical 
NMDAR is I(t) = Y^ =1 SigH(V(t)), where g is the single 
channel NMDAR conductance, and Si = 1 when recep- 
tor in open state when in closed and unbound states. 
To simplify the notation we can rewrite equation \J\ as 
H = V + ViQ{t-t post )e-^- t '"" t ^ TB , where V = a + bV R 
and V\ — bVs- 

Calcium concentration is assumed to be governed by 
equation ^ therefore the average calcium concentration 
will be: 



({Ca](t)) = ge-*'T f dt'e 1 '^ (j^ Si (t')\ H(t'), (32) 



where ( ) denotes averaging over the probability distri- 
bution of Si(t). 

However, (i>;(t)) = Z (sj(i)) = Zp°(t), thus the cal- 
cium concentration given Z NMDAR is: 



([Ca](t)) z = Zge- t/T / dt'e* ' l T p°(t')H{t'). (33) 



Expression is similar to equation [5] since I(t) = 
ZgH(t)p°(t). We have calculated this integral in section 
[n]for the different cases (pre, pre-post, post-pre). These 
assumptions are equivalent to the simple Markov model 
assumed here. This result is shown in equations 1321 and 
ESI We note that ([Ca](t)) z = Z {[Ca\(t)) r 

The variance for Z NMDAR, has the form a%(t) = 
([Ca] 2 (i)) z — {[Ca](t)) z . The second moment of the cal- 
cium concentration for Z NMDAR is: 



(29) ([Ca] 2 {t)) z =g 2 e- 2t l T df ' dt" V ' ' 7 e 4 ' ' / T H{t')H(t") 



i z z 
\i=l j=l 



2 e~ 2t/T / dt'dt"e t ' /T e t "f T H(t')H{t") 



-OO 

Z 

x t"), (34) 



where p°j(t, t') is the joint probability of having channel i 
open at time t and channel j open at time if . We assume 
independent identical channels, therefore for i ^ j we 
have p°°(t,t') = p°{t)p°(t'), and when i = j, p°°{t,t') is 
defined as p°°(t,t'). We therefore have that: 



([Ca] 2 {t)) z =g 2 e- 2t ' T / dt' dt" e*' ' T e*" ! T H(t')H{t") 

J — OO 

x [Z(Z-l)p°(t>)p°(t") + Zp°°(t\t")] 

=g 2 Z(Z-l)([Ca](t)) 2 +g 2 Ze- 2t ^ 

x f dt'dt"e t ' /T e t " /T H(t')H(t")p 00 (t',t") 

J — OO 

(35) 

Thus the variance of calcium concentration for Z NM- 
DAR is: 

a 2 (t) =g 2 Ze- 2t/T f dt 1 dt" / /r e l " l T H{t')H{t") 
x \p°°(t',t")-P°(t')p°(t")} 
=g 2 Ze- 2t / T f dt'dt"H{t')H{t")e t ' /T e t " /T p 00 (t',t") 



Z([Ca](t)) 2 . 



(36) 



We can use the Bayes rule to rewrite p°°(t,t') = 
p 00 (t\t')p°(t'). Thus there are two types of variables 
which we need to calculate from the Markov process, 
p°(t) and p°°(t\t'). The solution for p°(t) is given by 
equation [!n] We define the time t = as the beginning 
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(a) At = +10ms 
mean([Ca]) 



(b) At = -10ms 



100 200 300 
time (ms) 



(c) A t = +60 ms 




time (ms) 



100 200 
time (ms) 



FIG. 7: Mean calcium transients (solid line) and their standard deviation (dash-dot line) for different values of At, with Z = 10 
and /i = 0.5. (a) At = 10 ms, (b) At = —10 ms, and (c) At = 60 ms. The remaining parameters as in FIG. [^1 Note that the 
peaks of the calcium transients differ from that of FIG. |2| because fj, and Z are different. 



of the falling phase (G = 0), or t pre — 0. Further, we as- 
sume an instantaneous rise time. So for t < 0, p°{t) = 
and at t = p°(0) = fi. Thus for t > p°(t) = fie^ 1 /^ . 

Similarly, p°°(t\t') is the probability of being in open 
state at time t given that the channel was in open state 
at time t' . We calculate this for the case t, t' > 0. If t' > t 
this can be solved in a similar manner to p°, since in this 
regime it is a stationary Markov process, it must be a 
function of (t-f'), and at t = t, p°°\t = t'\t') = 1. During 
the falling phase there is one directional movement from 
o to u, since G = 0. Therefore, if we know that at time 
t' the channel is in an open state it must have been in an 
open state, at any time t < t' . Therefore: 



\t\t') = 



-(t-t')/TN 



t > t' 

t < a 



(37) 



Note that this holds for (t, t') > 0, and that by defini- 
tion at times smaller than zero the receptor is unbound. 
Hence: 



p°°(t,t') 



[i ■ e 
fi ■ e 



-t/TN 

-t'/rs 



t > t' 
t < i! 



(38) 



This expression can now be substituted into equation 
1361 and by using the linear approximation of TL we obtain 
an expression for erf (i) that can be computed analyti- 
cally. To aid in calculation we denote: 



a 2 z (t) = ZM 2 e- 2t / T (V* ai + 2V V ia2 + V 2 a 3 ) 



Z([Ca}(t)} 2 (39) 



where: 



x x = j dt' dt"p°°{t\t")e t '/ T e t " /T Q(t) (40a) 
Jo Jo 



OL2 



dt' 



*-(t'-t P ° St )/TB p OO 



(t', t") (40b) 



03=/ dt' dC'e-^-'^'/^e""-'"")/^/"^,*") 
Jo Jo 

t'/rj'/rr,,,, 



x e ' e 



0(t' - t post )Q(t" - t post ) 
(40c) 



In appendix ^ we analytically calculate the form of 
the terms oti, a.2 and Although the analytical form 
is complex, the consequences are simple and significant. 
As expected, CVz — CVi/V~Z. Thus the variability is 
significant only at relatively low Z. In FIG. we show 
([Ca]{t)) z and a z , for At = -10, 10 and 60 ms. In all 
cases both the mean and az change over time. The peak 
of ([Ca](t)) z for At = —10 ms and At = 60 ms are 
similar. However, the magnitude of the variability, a z , is 
significantly larger at At = 60 ms than at At = — 10 ms. 
To quantify the dependence of the relative variability as 
a function of At, we measure the coefficient of variation, 
CVz = ozl {[Ca](t)) z , where both [Ca](t) z and a z {t) 
are measured at the time t where ([Ca](t)) z is maximal. 

When At < the CV is low and decreases as At — > 
from below. For At > the CV increases as At increases 
(FIG |HJ) ■ It is interesting to compare the CV for values 
of At < and At > which have similar peak calcium 
levels. For example for At = —10 ms, CViq = 0.34 
and for At = 60 ms, which has a similar peak calcium 
level, CV W = 0.51. Therefore CVio(60)/CFio(-10) « 
1.5. Since CVz = CA\/\/~Z, for this set of parameters, 
but for every Z CV Z (60)/CV Z (-10) w 1.5. 

The variability of calcium transients also depends on fi 
(FIG.|5J), decreasing as /j, increases, for all At. The vari- 
able t-t is the probability of glutamate binding to post- 
synaptic receptors given a presynaptic spike, and can be 
taken to be the presynaptic probability of release. Thus, 
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FIG. 8: Coefficient of variation for Z = 10 as a function 
of At. Notice that the variability decreases with increasing 
values of r N . Shown are the plots for r N = 50 ms (circle), tn 
= 75 ms (square) and r N = 100 ms (diamond). 



CV 




0.2 0.4 0.6 0.8 
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FIG. 9: Coefficient of variation for Z = 10 as a function of 
fi, for three different values of At: At = — 10 ms (diamond), 
At = 10 ms (square) and At = 60 ms (circle) . As /i increases, 
CV decreases. 



as the presynaptic probability of release increases, CV 
decreases. 



VI. DISCUSSION 

Calcium transients due to NMDA currents are believed 
to play a major role in many forms of synaptic plasticity. 
We have recently shown that a model based on these 
transients can account for various different plasticity- 
induction protocols 28] . In this paper, we compute ana- 
lytically the calcium dynamics evoked by pairs of pre- and 
postsynaptic spikes under different conditions, so that 
the variables that control calcium transients and impact 
synaptic plasticity can be investigated. 

We showed that the peak of the calcium transients de- 
pends on the At: for At > 0, the peak concentration 
decays with the time constant t n of the NMDAR, and 
for At < 0, it decays with the time constant r B of the 
BPAP. Therefore, if spike time dependent plasticity, in- 



deed depends on the calcium transients, t n and r B will 
determine the width of the pre-post LTP and the post-pre 
LTD windows, respectively. These results were confirmed 
by simulations, in which the approximations we used in 
the calculations were relaxed. The At dependence of the 
peak calcium transient is similar in the more complex and 
realistic case where the BPAP dynamics is composed of 
a fast and a slow component: the difference lies in the 
At < case, which now depends on the combination 
of the two time constants. The two-component BPAP 
allows a sharper transition between the post-before-pre 
and the pre-before-post STDP window, which is further 
enhanced by a fast calcium dynamics. We also showed 
that, besides the timing of the pre and post-spikes, the 
peak values of the calcium transients also depend on the 
frequency v of the pre and postsynaptic conditioning, in- 
creasing for greater values of v. However, at higher fre- 
quencies, the dependence on At decreases. The amount 
of calcium build-up at a given frequency depends on the 
system parameters; slower time constants results in more 
temporal integration. Temporal integration of calcium 
transients, as described here, can explain why, in some 
cases, the induction of STDP is frequency-dependent. 

One of the predictions of the unified calcium model is 
that there exists a pre-before-post LTD for values of At 
greater than the ones that elicit LTP [23 • This has been 
shown experimentally by some investigators |24j : others 
have only placed a small number of data points in this 
region (a S3- I n au published STDP experiments, there 
is a large variability in the magnitude and the sign of 
synaptic plasticity across the different values of At. It 
would therefore be difficult to assess the existence of this 
type of LTD without a large amount of data. The vari- 
ability encountered in the experiments indicates that it 
should be important to examine the stochastic properties 
of the calcium transients, in addition to their mean. 

We calculated the variance of the calcium transients 
and showed that for a small number Z of NMDAR, it 
can be quite significant. There are indications that the 
number of NMDAR in the dendrite is indeed small (« 10) 
|25j . Further, we showed that the CV increases monoton- 
ically with At, for At > 0. These fluctuations can have 
significant implications on the outcome of downstream 
processes dependent on calcium, such as phosphorylation 
or other second-messenger cascades upon which plasticity 
relies. 

Previously, several authors have simulated and ana- 
lyzed dynamical models of synaptic plasticity that can 
lead to STDP [1 El El lH Hj. These models typi- 
cally bypass detailed descriptions of the physiology and 
biochemistry underlying synaptic plasticity. Recently, 
bio phy sical models have been receiving increasing focus 
HI LLil liii HE • These more mechanistic approaches 
have allowed one to show that various different induc- 
tion protocols can be accounted for by a single set of 
assumptions regarding the dependence of plasticity on 
calcium concentrations [2Sj. In these models, plasticity is 
a highly non-linear function of calcium. Therefore, flue- 
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tuations of calcium transients will not only add spread 
to the resulting LTD and LTP, but will also shift their 
mean (the mean of the solution is not equivalent to the 
solution given the mean calcium level). For example, in 
the unified calcium model |28j . the sign and magnitude 
of synaptic plasticity is determined by a saturated U- 
shaped function of calcium concentration. When calcium 
falls below a threshold 9d, no plasticity occurs, when it 
falls between 9d and 9 pi 9 p > Od, LTD is induced, and for 
calcium above 9 p , LTP is induced. Thus, in the LTD re- 
gion, if the spread of the calcium is of the order of 9 p — 0d, 
occasional LTP can still happen. We showed an explicit 
example for the cases of At — — 10 ms and At = 60 ms, 
where the mean peak calcium concentrations are similar, 
but the variability for the latter is larger. If this mean 
calcium is in the LTD range, it is likely that the observed 
phenomenon is that the LTD is stronger and more robust 
for At = — 10 ms than for At = 60 ms. 

This work is based on several simplifying assumptions. 
This allows us a qualitative understanding of the aspects 
of calcium dynamics that are crucial for the induction 
of synaptic plasticity. Incorporation of more realistic 
features, such as those assumed in the unified calcium 
model, would probably yield quantitatively more precise 
results. For example, a double-component NMDAR de- 
cay would allow for more calcium build-up without sac- 
rificing the temporal asymmetry between pre-post and 
post-pre conditions. However, we believe that, despite 
the constraints of analytical tractability, this study can 
serve as a basis for future models that aim to formalize 
the biophysical basis of synaptic plasticity. 



' dt' f dt"e- t "'^e t '' T e t "' T 
o Jt' 

+ f dt' f dt"e- t '^e t '/ T e t "/ T 
Jo Jo 



Q(t) 



+ [d^r(e^-l) d , 



©(*) 



3*/ T2 t 2 t ((e*/ T - l) -T 2 ( dt'e*'/ 



+ t dt'e 1 '!^ -t dt'e*'' 1 
Jo Jo 

e t/Tl (t 2 T - T 2 T4 + TT A ) 



e(t) 



-2e t/T2 r 2 r + (t 2 t 4 - rr 4 + t 2 t)1 Q(t), (Al) 



where t 4 1 = 2r 1 — r N x . 



a 2 =e 



t post /r B 



T 5 (T-T 2 )(e^-e f > 5 )- 



-TT 3 (e 



t/T3 _ J/T 3 \ 



T 2 T 6 e* /T2 (e t/T6 -e i/Te ) 



x G(t - 1) 



(A2) 



where T3 is as defined before, r 5 1 = 2r 1 — r N 1 — r B x , 



j-post 





t pOSt > Q 
t pOSt < Q 



(A3) 



Finally, 
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APPENDIX B: NUMERICAL ESTIMATION OF 
THE ERROR RESULTING FROM THE LINEAR 
APPROXIMATION OF H 

In the simulations, we used the following form for the 
voltage-dependence of the NMDAR.: 



APPENDIX A: VARIANCE CALCULATION 

Using the definition in equation I40al through I40cl we 
obtain: 



n(v) - e _ Z v 2V (bi) 

where V rev — 130 mV, fit performed in the interval 
[-70,-10] mV. 

Such functional dependence captures the qualitative 
non-linear dependence of the NMDAR on the voltage, al- 
though the reversal potential for calcium channels might 
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FIG. 10: The error estimation for the linear approximation of 
the 7i-function. a) The original 7i-function (dotted line) and 
its linear approximation (solid line) in the interval [—70, —10] 
mV. b), c) and d) The calcium transients elicited by the pre- 
only, post-pre and pre-post conditions, respectively, for the 
full expression (dotted line) and for the linear approximation 
(solid line). 



be inaccurate at the values close to it. In our analysis, 
we approximated H. with a linear dependence on volt- 
age. The fit used is shown in FIG. HOfa . We see that 
over the relevant range ([-70 0]) this is a relatively good 
approximation. 



In panels b-d of FIG. EH we compare the numerically 
extracted calcium transients (dashed line), using the non- 
linear 7i and the analytical results obtained for the linear 
approximation (solid line). At rest (FIG. 110b). the lin- 
ear approximation significantly underestimates the mag- 
nitudes of the transients, because this is where the linear 
fit departs the most from the true curve. However, in the 
working regime ([-10, 10] mV), the approximation can be 
considered adequate. 



Let equation IB1I be rewritten as Ti.(v) = a + bV as 
in equation [7| The parameters of linearization are a = 
0.1031/g and b = 0.0015/0, where g = -10~ 3 is the 
average NMDA conductance. 
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